!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine regrid (current)
!==================================================================================================================================|

   use parameters
   use block_define
   use interface_block_generate
!   use interface_block_destroy
!   use interface_tvd_mmc_x
!   use interface_tvd_mmc_y
   use interface_hilbert_space_filling_2d
   implicit none

   type(block), pointer, intent(inout) :: current

   integer(4) :: i, i_block, j, j_block, i_neigh, j_neigh, ip1, im1, jp1, jm1
   integer(4) :: flag, flag_destroy, n_blocks_tmp, n_gene_tmp
   integer(4), dimension(:, :, :), allocatable :: gene_tmp

   real(8) :: lv_exp, num_tmp, thral1
!   real(8), dimension(n_neighs, n_dims) :: len_p_bnd
   real(8), dimension(n_thr) :: dva, dva2, d2va2
   real(8), dimension(nx_block, ny_block) :: vx, vy, vz, pr, te
!   real(8), dimension(nx_block, ny_block) :: ro_c, mx_c, my_c, mz_c, bx_c, by_c, bz_c, en_c, gx_c, gy_c
!   real(8), dimension(nx_block, ny_block) :: ro_p, mx_p, my_p, mz_p, bx_p, by_p, bz_p, en_p, gx_p, gy_p
   type(block), pointer :: child, neigh, next, framework_next, pointer_tmp

!----------------------------------------------------------------------------------------------------------------------------------|
!  error estimation (flag = 1 then block generation, flag = -1 then block destruction and flag = 0 for nothing)
!  Richardson extrapolation: regrid_flag = 1
!  Ziegler criterion: regrid_flag = 2
!----------------------------------------------------------------------------------------------------------------------------------|
   flag = 0

   if (current%lv .lt. n_levels) then

      flag = -1

!      if (regrid_flag .eq. 1) then
!         if (current%lv .eq. 1) then
!            loop1: do j = ng + 1, current%ny - ng
!            do i = ng + 1, current%nx - ng
!               ip1 = i + 1
!               im1 = i - 1
!               jp1 = j + 1
!               jm1 = j - 1
!               dva(1) = (abs(current%ro(ip1, j) - current%ro(im1, j)) +                                                           &
!                  abs(current%ro(i, jp1) - current%ro(i, jm1))) / (abs(current%ro(i, j)) + eps)
!               dva(2) = (abs(current%mx(ip1, j) - current%mx(im1, j)) +                                                           &
!                  abs(current%mx(i, jp1) - current%mx(i, jm1))) / (abs(current%mx(i, j)) + eps)
!               dva(3) = (abs(current%my(ip1, j) - current%my(im1, j)) +                                                           &
!                  abs(current%my(i, jp1) - current%my(i, jm1))) / (abs(current%my(i, j)) + eps)
!               dva(4) = (abs(current%mz(ip1, j) - current%mz(im1, j)) +                                                           &
!                  abs(current%mz(i, jp1) - current%mz(i, jm1))) / (abs(current%mz(i, j)) + eps)
!               dva(5) = (abs(current%bx(ip1, j) - current%bx(im1, j)) +                                                           &
!                  abs(current%bx(i, jp1) - current%bx(i, jm1))) / (abs(current%bx(i, j)) + eps)
!               dva(6) = (abs(current%by(ip1, j) - current%by(im1, j)) +                                                           &
!                  abs(current%by(i, jp1) - current%by(i, jm1))) / (abs(current%by(i, j)) + eps)
!               dva(7) = (abs(current%bz(ip1, j) - current%bz(im1, j)) +                                                           &
!                  abs(current%bz(i, jp1) - current%bz(i, jm1))) / (abs(current%bz(i, j)) + eps)
!               dva(8) = (abs(current%en(ip1, j) - current%en(im1, j)) +                                                           &
!                  abs(current%en(i, jp1) - current%en(i, jm1))) / (abs(current%en(i, j)) + eps)
!               if (any (dva .gt. thr_va0)) then
!                  flag = 1
!                  len_p_bnd(1, 1) = current%x(i) - current%x(1)
!                  len_p_bnd(2, 1) = current%x(current%nx) - current%x(i)
!                  len_p_bnd(1, 2) = current%y(j) - current%y(1)
!                  len_p_bnd(2, 2) = current%y(current%ny) - current%y(j)
!                  len_regrid = min (len_regrid, minval (len_p_bnd))
!                  exit loop1
!               else if (any (dva .gt. thr_va0 * thr_de)) then
!                  flag = 0
!               endif
!            enddo
!            enddo loop1
!         else
!            parent => current%parent
!            do j = 1, parent%ny
!            do i = 1, parent%nx
!               ro_p(i, j) = parent%ro(i, j)
!               mx_p(i, j) = parent%mx(i, j)
!               my_p(i, j) = parent%my(i, j)
!               mz_p(i, j) = parent%mz(i, j)
!               bx_p(i, j) = parent%bx(i, j)
!               by_p(i, j) = parent%by(i, j)
!               bz_p(i, j) = parent%bz(i, j)
!               en_p(i, j) = parent%en(i, j)
!               gx_p(i, j) = parent%gx(i, j)
!               gy_p(i, j) = parent%gy(i, j)
!            enddo
!            enddo

!            call tvd_mmc_x (ro_p, mx_p, my_p, mz_p, bx_p, by_p, bz_p, en_p, gx_p,                                                  &
!               parent%nx, parent%ny, ng, n_neighs, parent%dx, parent%dt, max_speed, gm, div_clean_flag)
!            call tvd_mmc_y (ro_p, mx_p, my_p, mz_p, bx_p, by_p, bz_p, en_p, gy_p,                                                  &
!               parent%nx, parent%ny, ng, n_neighs, parent%dy, parent%dt, max_speed, gm, div_clean_flag)

!            do j = 1, current%ny
!            do i = 1, current%nx
!               ro_c(i, j) = current%ro(i, j)
!               mx_c(i, j) = current%mx(i, j)
!               my_c(i, j) = current%my(i, j)
!               mz_c(i, j) = current%mz(i, j)
!               bx_c(i, j) = current%bx(i, j)
!               by_c(i, j) = current%by(i, j)
!               bz_c(i, j) = current%bz(i, j)
!               en_c(i, j) = current%en(i, j)
!               gx_c(i, j) = current%gx(i, j)
!               gy_c(i, j) = current%gy(i, j)
!            enddo
!            enddo

!            call tvd_mmc_x (ro_c, mx_c, my_c, mz_c, bx_c, by_c, bz_c, en_c, gx_p,                                                  &
!               current%nx, current%ny, ng, n_neighs, current%dx, current%dt, max_speed, gm, div_clean_flag)
!            call tvd_mmc_y (ro_c, mx_c, my_c, mz_c, bx_c, by_c, bz_c, en_c, gy_p,                                                  &
!               current%nx, current%ny, ng, n_neighs, current%dy, current%dt, max_speed, gm, div_clean_flag)

!            lv_exp = 2 ** (thr_ex * parent%lv)
!            nx_h = nx_no_bnd / parent%nc_x
!            ny_h = ny_no_bnd / parent%nc_y
!            loop2: do j = 1, ny_h
!            do i = 1, nx_h
!               i_p = i + nx_h * (current%p_pos(1) - 1) + ng
!               i_c = i * 2 + ng
!               icp1 = i_c - 1
!               j_p = j + ny_h * (current%p_pos(2) - 1) + ng
!               j_c = j * 2 + ng
!               jcp1 = j_c - 1
!               dva(1) = abs(ro_p(i_p, j_p) - (ro_c(icp1, j_c) + ro_c(i_c, j_c) + ro_c(icp1, jcp1) + ro_c(i_c, jcp1)) / 4.0d0) /   &
!                  (abs(ro_p(i_p, j_p) + eps)) * lv_exp
!               dva(2) = abs(mx_p(i_p, j_p) - (mx_c(icp1, j_c) + mx_c(i_c, j_c) + mx_c(icp1, jcp1) + mx_c(i_c, jcp1)) / 4.0d0) /   &
!                  (abs(mx_p(i_p, j_p) + eps)) * lv_exp
!               dva(3) = abs(my_p(i_p, j_p) - (my_c(icp1, j_c) + my_c(i_c, j_c) + my_c(icp1, jcp1) + my_c(i_c, jcp1)) / 4.0d0) /   &
!                  (abs(my_p(i_p, j_p) + eps)) * lv_exp
!               dva(4) = abs(mz_p(i_p, j_p) - (mz_c(icp1, j_c) + mz_c(i_c, j_c) + mz_c(icp1, jcp1) + mz_c(i_c, jcp1)) / 4.0d0) /   &
!                  (abs(mz_p(i_p, j_p) + eps)) * lv_exp
!               dva(5) = abs(bx_p(i_p, j_p) - (bx_c(icp1, j_c) + bx_c(i_c, j_c) + bx_c(icp1, jcp1) + bx_c(i_c, jcp1)) / 4.0d0) /   &
!                  (abs(bx_p(i_p, j_p) + eps)) * lv_exp
!               dva(6) = abs(by_p(i_p, j_p) - (by_c(icp1, j_c) + by_c(i_c, j_c) + by_c(icp1, jcp1) + by_c(i_c, jcp1)) / 4.0d0) /   &
!                  (abs(by_p(i_p, j_p) + eps)) * lv_exp
!               dva(7) = abs(bz_p(i_p, j_p) - (bz_c(icp1, j_c) + bz_c(i_c, j_c) + bz_c(icp1, jcp1) + bz_c(i_c, jcp1)) / 4.0d0) /   &
!                  (abs(bz_p(i_p, j_p) + eps)) * lv_exp
!               dva(8) = abs(en_p(i_p, j_p) - (en_c(icp1, j_c) + en_c(i_c, j_c) + en_c(icp1, jcp1) + en_c(i_c, jcp1)) / 4.0d0) /   &
!                  (abs(en_p(i_p, j_p) + eps)) * lv_exp
!               if (any (dva .gt. thr_va)) then
!                  flag = 1
!                  len_p_bnd(1, 1) = parent%x(i_p) - parent%x(1)
!                  len_p_bnd(2, 1) = parent%x(parent%nx) - parent%x(i_p)
!                  len_p_bnd(1, 2) = parent%y(j_p) - parent%y(1)
!                  len_p_bnd(2, 2) = parent%y(parent%ny) - parent%y(j_p)
!                  len_regrid = min (len_regrid, minval (len_p_bnd))
!                  exit loop2
!               else if (any (dva .gt. thr_va * thr_de)) then
!                  flag = 0
!               endif
!            enddo
!            enddo loop2
!         endif
!      elseif (regrid_flag .eq. 2) then
      lv_exp = dble(ratio) ** (thr_ex * (current%lv - 1))
      num_tmp = sqrt(2.0d0) / 2.0d0
      thral1 = 1.0d0 - thr_al
      call variables_transform (vx, vy, vz, pr, te,                                                                                &
         current%ro, current%mx, current%my, current%mz, current%bx, current%by, current%bz, current%en, gm, nx_block, ny_block, 1)
      loop3: do j = 2, ny_block - 1
         jp1 = j + 1
         jm1 = j - 1
         do i = 2, nx_block - 1
            ip1 = i + 1
            im1 = i - 1

            dva2(1) = num_tmp * sqrt((current%ro(ip1, j) - current%ro(im1, j)) ** 2 +                                              &
                                      (current%ro(i, jp1) - current%ro(i, jm1)) ** 2)
            dva2(2) = num_tmp * sqrt((vx(ip1, j) - vx(im1, j)) ** 2 + (vx(i, jp1) - vx(i, jm1)) ** 2)
            dva2(3) = num_tmp * sqrt((vy(ip1, j) - vy(im1, j)) ** 2 + (vy(i, jp1) - vy(i, jm1)) ** 2)
            dva2(4) = num_tmp * sqrt((vz(ip1, j) - vz(im1, j)) ** 2 + (vz(i, jp1) - vz(i, jm1)) ** 2)
            dva2(5) = num_tmp * sqrt((current%bx(ip1, j) - current%bx(im1, j)) ** 2 +                                              &
                                      (current%bx(i, jp1) - current%bx(i, jm1)) ** 2)
            dva2(6) = num_tmp * sqrt((current%by(ip1, j) - current%by(im1, j)) ** 2 +                                              &
                                      (current%by(i, jp1) - current%by(i, jm1)) ** 2)
            dva2(7) = num_tmp * sqrt((current%bz(ip1, j) - current%bz(im1, j)) ** 2 +                                              &
                                      (current%bz(i, jp1) - current%bz(i, jm1)) ** 2)
            dva2(8) = num_tmp * sqrt((pr(ip1, j) - pr(im1, j)) ** 2 + (pr(i, jp1) - pr(i, jm1)) ** 2)
            dva2(9) = num_tmp * sqrt((te(ip1, j) - te(im1, j)) ** 2 + (te(i, jp1) - te(i, jm1)) ** 2)

            d2va2(1) = sqrt((current%ro(ip1, j) - 2.0d0 * current%ro(i, j) + current%ro(im1, j)) ** 2 +                            &
                             (current%ro(i, jp1) - 2.0d0 * current%ro(i, j) + current%ro(i, jm1)) ** 2)
            d2va2(2) = sqrt((vx(ip1, j) - 2.0d0 * vx(i, j) + vx(im1, j)) ** 2 + (vx(i, jp1) - 2.0d0 * vx(i, j) + vx(i, jm1)) ** 2)
            d2va2(3) = sqrt((vy(ip1, j) - 2.0d0 * vy(i, j) + vy(im1, j)) ** 2 + (vy(i, jp1) - 2.0d0 * vy(i, j) + vy(i, jm1)) ** 2)
            d2va2(4) = sqrt((vz(ip1, j) - 2.0d0 * vz(i, j) + vz(im1, j)) ** 2 + (vz(i, jp1) - 2.0d0 * vz(i, j) + vz(i, jm1)) ** 2)
            d2va2(5) = sqrt((current%bx(ip1, j) - 2.0d0 * current%bx(i, j) + current%bx(im1, j)) ** 2 +                            &
                             (current%bx(i, jp1) - 2.0d0 * current%bx(i, j) + current%bx(i, jm1)) ** 2)
            d2va2(6) = sqrt((current%by(ip1, j) - 2.0d0 * current%by(i, j) + current%by(im1, j)) ** 2 +                            &
                             (current%by(i, jp1) - 2.0d0 * current%by(i, j) + current%by(i, jm1)) ** 2)
            d2va2(7) = sqrt((current%bz(ip1, j) - 2.0d0 * current%bz(i, j) + current%bz(im1, j)) ** 2 +                            &
                             (current%bz(i, jp1) - 2.0d0 * current%bz(i, j) + current%bz(i, jm1)) ** 2)
            d2va2(8) = sqrt((pr(ip1, j) - 2.0d0 * pr(i, j) + pr(im1, j)) ** 2 + (pr(i, jp1) - 2.0d0 * pr(i, j) + pr(i, jm1)) ** 2)
            d2va2(9) = sqrt((te(ip1, j) - 2.0d0 * te(i, j) + te(im1, j)) ** 2 + (te(i, jp1) - 2.0d0 * te(i, j) + te(i, jm1)) ** 2)

            dva(1) = (thr_al * dva2(1) / (abs(current%ro(i, j)) + eps) +                                                           &
                      thral1 * d2va2(1) / (dva2(1) + thr_fl * (abs(current%ro(i, j)) + eps))) * lv_exp
            dva(2) = (thr_al * dva2(2) / (abs(vx(i, j)) + 0.1d0) +                                                                 &
                      thral1 * d2va2(2) / (dva2(2) + thr_fl * (abs(vx(i, j)) + 0.1d0))) * lv_exp
            dva(3) = (thr_al * dva2(3) / (abs(vy(i, j)) + 0.1d0) +                                                                 &
                      thral1 * d2va2(3) / (dva2(3) + thr_fl * (abs(vy(i, j)) + 0.1d0))) * lv_exp
            dva(4) = (thr_al * dva2(4) / (abs(vz(i, j)) + 0.1d0) +                                                                 &
                      thral1 * d2va2(4) / (dva2(4) + thr_fl * (abs(vz(i, j)) + 0.1d0))) * lv_exp
            dva(5) = (thr_al * dva2(5) / (abs(current%bx(i, j)) + 0.1d0) +                                                         &
                      thral1 * d2va2(5) / (dva2(5) + thr_fl * (abs(current%bx(i, j)) + 0.1d0))) * lv_exp
            dva(6) = (thr_al * dva2(6) / (abs(current%by(i, j)) + 0.1d0) +                                                         &
                      thral1 * d2va2(6) / (dva2(6) + thr_fl * (abs(current%by(i, j)) + 0.1d0))) * lv_exp
            dva(7) = (thr_al * dva2(7) / (abs(current%bz(i, j)) + 0.1d0) +                                                         &
                      thral1 * d2va2(7) / (dva2(7) + thr_fl * (abs(current%bz(i, j)) + 0.1d0))) * lv_exp
            dva(8) = (thr_al * dva2(8) / (abs(pr(i, j)) + eps) +                                                                   &
               thral1 * d2va2(8) / (dva2(8) + thr_fl * (abs(pr(i, j)) + eps))) * lv_exp
            dva(9) = (thr_al * dva2(9) / (abs(te(i, j)) + eps) +                                                                   &
               thral1 * d2va2(9) / (dva2(9) + thr_fl * (abs(te(i, j)) + eps))) * lv_exp

            if (any (dva .gt. thr_va0)) then
               flag = 1
               exit loop3
            else if (any (dva .gt. thr_va0 * thr_de)) then
               flag = 0
            endif
         enddo
      enddo loop3
!      endif
   endif

!----------------------------------------------------------------------------------------------------------------------------------|
!  grid generation or destruction
!----------------------------------------------------------------------------------------------------------------------------------|
   n_blocks_tmp = n_blocks(cpu_id_p1)
   if (flag .eq. 1) then

      if (current%nc .eq. 0 .and. current%lv .lt. n_levels) then
         i_gene(cpu_id_p1) = i_gene(cpu_id_p1) + 1
!----------------------------------------------------------------------------------------------------------------------------------|
!        for the case i_gene(cpu_id_p1) is greater than the n_gene
!----------------------------------------------------------------------------------------------------------------------------------|
         if (i_gene(cpu_id_p1) .gt. n_gene) then
            allocate(gene_tmp(n_dims, n_levels, n_gene))
            gene_tmp = gene
            n_gene_tmp = n_gene
            n_gene = n_gene + ini_block_x * ini_block_y
            deallocate(gene)
            allocate(gene(n_dims, n_levels, n_gene))
            gene(:, :, 1 : n_gene_tmp) = gene_tmp
            deallocate(gene_tmp)
         endif
         gene(:, :, i_gene(cpu_id_p1)) = current%id
         call block_generate (current, current%cpu_id)
         next => current%next
         framework_next => current%framework_next
         pointer_tmp => current
         call hilbert_space_filling_2d (current, pointer_tmp, current%start_pos, current%end_pos, current%start_point,             &
            ratio_order, current%cpu_id)
         next%last => pointer_tmp
         pointer_tmp%next => next
         framework_next%framework_last => pointer_tmp
         pointer_tmp%framework_next => framework_next
         max_level = max(max_level, current%lv + 1)
      endif

   elseif (flag .eq. -1) then

      flag_destroy = 1
      do j_block = 1, current%nc_y
      do i_block = 1, current%nc_x
         child => current%child(i_block, j_block)%pblock
         do j_neigh = -1, 1
         do i_neigh = -1, 1
            neigh => child%neigh(i_neigh, j_neigh)%pblock
            if (associated (neigh)) then
               if (all (neigh%parent%id .eq. child%parent%id)) cycle
               if (neigh%lv .eq. child%lv .and. neigh%nc .ne. 0) flag_destroy = 0
            endif
         enddo
         enddo
      enddo
      enddo

      if (flag_destroy .eq. 1 .and. current%nc .ne. 0) then
         i_gene(cpu_id_p1) = i_gene(cpu_id_p1) + 1
!----------------------------------------------------------------------------------------------------------------------------------|
!        for the case i_gene(cpu_id_p1) is greater than the n_gene
!----------------------------------------------------------------------------------------------------------------------------------|
         if (i_gene(cpu_id_p1) .gt. n_gene) then
            allocate(gene_tmp(n_dims, n_levels, n_gene))
            gene_tmp = gene
            n_gene_tmp = n_gene
            n_gene = n_gene + ini_block_x * ini_block_y
            deallocate(gene)
            allocate(gene(n_dims, n_levels, n_gene))
            gene(:, :, 1 : n_gene_tmp) = gene_tmp
            deallocate(gene_tmp)
         endif
         gene(:, :, i_gene(cpu_id_p1)) = - current%id
      endif

   endif

!----------------------------------------------------------------------------------------------------------------------------------|
!  reset n_blocks
!----------------------------------------------------------------------------------------------------------------------------------|
   n_blocks(cpu_id_p1) = n_blocks_tmp

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine regrid
